
#include "errlog.h"
#include "grid3d.h"
#include "vecs.h"
#include "quadr.h"
#include "fe_bubble.h"

#include <assert.h>
#include <stdlib.h>
#include <memory.h>

#define BUB_EXTEND_ID           0x3FF11624      /* a random number
* large enough so it is likely
* to be unique
*/

/*
 * Initialize the bubble element
 **********************************
 * This routine augments the nodes array with the centroids of the
 * tetrahedra.  It also creates an extended connectivity table by adding
 * one node (the centroid) to each element.
 * Nothing in the connectivity table changes since the indices of the nodes at
 * the vertices of the tetrahedra do not change.  We only add more nodes at
 * the end.
 */
void InitBubble(grid3d_t *grid)
{
    /* don't use two bubble grids at the same time!! */
    static bub_extend_t bub_extend;
    elem_t *T;
    int el, j;

    assert(grid!=NULL);
    bub_extend.id=BUB_EXTEND_ID;
    /* allocate extended grid data */
    bub_extend.b_el=malloc(grid->NEL*sizeof(bub_elem_t));
    if(bub_extend.b_el==NULL)
        error_msg("Memory allocation error in InitBubble.\n");
    grid->extended= &bub_extend;
    /* allocate the new nodes */
    bub_extend.bubNN=grid->NN+grid->NEL;
#if FALSE
    nds=malloc(bub_extend.bubNN*sizeof(node_t));
    if(nds==NULL)
        error_msg("Memory allocation error in InitBubble.\n");
    memcpy(nds, grid->nodes, grid->NN*sizeof(node_t));
    free(grid->nodes);
    grid->nodes=nds;
#endif
    /* fill in the blanks */
    for(el=0; el<grid->NEL; el++)
    {
        T=grid->elems+el;
#if FALSE
        int i;
        node_t *nds;
        nds=grid->nodes+grid->NN+el;
        for(i=0; i<3; i++)
            nds->coords[i]=0.0;
        for(i=0; i<3; i++)
            for(j=0; j<4; j++)
                nds->coords[i]+=grid->nodes[T->nodes[j]].coords[i];
        for(i=0; i<3; i++)
            nds->coords[i]=nds->coords[i]*0.25;
        nds->bcond=BC_NON; /* not on boundary */
        nds->on_bdry=0;
#endif
        for(j=0; j<4; j++)
        {
            bub_extend.b_el[el].nodes[j]=T->nodes[j];
            bub_extend.b_el[el].dofs[j]=T->nodes[j];
            bub_extend.b_el[el].dofs[j+5]=T->nodes[j]+bub_extend.bubNN;
            bub_extend.b_el[el].dofs[j+10]=T->nodes[j]+bub_extend.bubNN*2;
        }
        bub_extend.b_el[el].nodes[4]=grid->NN+el;
        bub_extend.b_el[el].dofs[4]=grid->NN+el;
        bub_extend.b_el[el].dofs[4+5]=grid->NN+el+bub_extend.bubNN;
        bub_extend.b_el[el].dofs[4+10]=grid->NN+el+bub_extend.bubNN*2;
    }
}

static bub_extend_t *get_be(const grid3d_t *grid)
{
    bub_extend_t *be;
    assert(grid!=NULL);
    be=grid->extended;
    if((be==NULL)||(be->id!=BUB_EXTEND_ID))
        error_msg("fe_bubble.c: grid->extended is not a valid "
                  "bub_extend_t structure .\n");
    return be;
}

/* Finalize the bubble element
 *************************************
 * This routine simply releases the memory taken by the extended grid data
 */
void DoneBubble(grid3d_t *grid)
{
    bub_extend_t *bub_ext=get_be(grid);
    free(bub_ext->b_el);
}

int ndof_bub(const grid3d_t *grid)
{
    bub_extend_t *be=get_be(grid);
    return be->bubNN;
}

int ndof_bub_vec(const grid3d_t *grid)
{
    bub_extend_t *be=get_be(grid);
    return 3*be->bubNN;
}

/****************************************************************************/

void dof_coords_bub(const grid3d_t *grid, node_t *nodes,
                    int dof, double **coords)
{
    bub_extend_t *be=get_be(grid);
    if(dof>=be->bubNN)
        error_msg("%s: dof out of range.\n", __func__);
    if(coords)
    {
        if(dof<grid->NN)
            *coords=nodes[dof].coords;
        else
        {
            static double pt[3];
            ComputeCentroid(nodes, grid->elems[dof-grid->NN].nodes, pt);
            *coords=pt;
        }
    }
}

void dof_coords_bub_vec(const grid3d_t *grid, node_t *nodes,
                        int dof, double **coords)
{
    bub_extend_t *be=get_be(grid);
    if(dof>=be->bubNN*3)
        error_msg("%s: dof out of range.\n", __func__);
    dof_coords_bub(grid, nodes, dof%be->bubNN, coords);
}

/****************************************************************************/

unsigned dof_bcond_bub(const grid3d_t *grid, int dof)
{
    bub_extend_t *be=get_be(grid);
    if(dof>=be->bubNN)
        error_msg("%s: dof out of range.\n", __func__);
    if(dof<grid->NN)
        return grid->nodes[dof].bcond;
    else
        return BC_NON;
}

unsigned dof_bcond_bub_vec(const grid3d_t *grid, int dof)
{
    bub_extend_t *be=get_be(grid);
    if(dof>=3*be->bubNN)
        error_msg("%s: dof out of range.\n", __func__);
    return dof_bcond_bub(grid, dof%be->bubNN);
}

/****************************************************************************/

unsigned dof_bdry_bub(const grid3d_t *grid, int dof)
{
    bub_extend_t *be=get_be(grid);
    if(dof>=be->bubNN)
        error_msg("%s: dof out of range.\n", __func__);
    if(dof<grid->NN)
        return grid->nodes[dof].on_bdry;
    else
        return 0;
}

unsigned dof_bdry_bub_vec(const grid3d_t *grid, int dof)
{
    bub_extend_t *be=get_be(grid);
    if(dof>=be->bubNN*3)
        error_msg("%s: dof out of range.\n", __func__);
    return dof_bdry_bub(grid, dof%be->bubNN);
}

/****************************************************************************/

int nodes_bub(const grid3d_t *grid, int el, int **nodes)
{
    bub_extend_t *be=get_be(grid);
    if(nodes)
        *nodes=be->b_el[el].nodes;
    return 5;
}

int nodes_bub_vec(const grid3d_t *grid, int el, int **nodes)
{
    bub_extend_t *be=get_be(grid);
    if(nodes)
        *nodes=be->b_el[el].nodes;
    return 5;
}

/****************************************************************************/

int dofs_bub(const grid3d_t *grid, int el, int **dofs)
{
    bub_extend_t *be=get_be(grid);
    if(dofs)
        *dofs=be->b_el[el].dofs;
    return 5;
}

int dofs_bub_vec(const grid3d_t *grid, int el, int **dofs)
{
    bub_extend_t *be=get_be(grid);
    if(dofs)
        *dofs=be->b_el[el].dofs;
    return 15;
}




/****************************************************************************/
/*       local P1+bubble matrices */
/****************************************************************************/

/* define some factorial constants (n!=n*(n-1)*...*3*2*1) */
#define FAC2    2           /* 2! */
#define FAC3    6           /* 3! */
#define FAC5    120         /* 5! */
#define FAC7    5040        /* 7! */
#define FAC8    40320       /* 8! */
#define FAC9    362880      /* 9! */
#define FAC10   3628800     /* 10! */
#define FAC11   39916800    /* 11! */

void local_mass_bub(const grid3d_t *grid, const elem_t *T, double ***M)
{
    int i, j;
    static double A[25], *AA[5]={A, A+5, A+10, A+15, A+20};
    double TV;
    const double
    SIJ=FAC3/(double) FAC5,
            SII=FAC3*FAC2/(double) FAC5,
            SI4=(1<<8)*FAC2*FAC3/(double) FAC8,
            S44=(1<<20)*FAC3/(double) FAC11;
    if(M==NULL)
        return;
    else
        *M=AA;
    TV=ComputeVolume(grid->nodes, T->nodes);
    for(i=0; i<4; i++)
        for(j=0; j<4; j++)
            (*M)[i][j]=TV*(i==j?SII:SIJ);
    for(i=0; i<4; i++)
        (*M)[i][4]=(*M)[4][i]=TV*SI4;
    (*M)[4][4]=TV*S44;
}

void local_mass_bub_vec(const grid3d_t *grid, const elem_t *T, double ***M)
{
    int i, j;
    static double A[15*15],
            *AA[15]={
        A, A+15, A+30, A+45, A+60,
        A+75, A+90, A+105, A+120, A+135,
        A+150, A+165, A+180, A+195, A+210
    };
    double **lM;
    if(M==NULL)
        return;
    else
        *M=AA;
    memset(A, 0, sizeof(A));
    local_mass_bub(grid, T, &lM);
    for(i=0; i<5; i++) for(j=0; j<5; j++)
            (*M)[i][j]=(*M)[i+5][j+5]=(*M)[i+10][j+10]=lM[i][j];
    return;
}

static void tensor_grad_bub_grad_bub(const grid3d_t *grid, const elem_t *T,
                                     double *****_G);

void local_curl_bub_curl_bub(const grid3d_t *grid, const elem_t *T, double ***R)
/*
 *  This is the matrix  with entries
 *     R[i][j] = \int (curl b_i).(curl b_j),   0<= i,j < 15,
 *  where
 *     b_i = [phi_i, 0, 0]^T, for  0<= i <5
 *     b_i = [0, phi_i, 0]^T, for  5<= i <10
 *     b_i = [0, 0, phi_i]^T, for 10<= i <15
 */
{
    int i, j;
    double ****G;

    static double B[225], *BB[15]={
        B, B+15, B+30, B+45, B+60,
        B+75, B+90, B+105, B+120, B+135,
        B+150, B+165, B+180, B+195, B+210
    };
    if(R==NULL)
        return;
    else
        *R=BB;
    tensor_grad_bub_grad_bub(grid, T, &G);
    memset(B, 0, sizeof(B));
    for(i=0; i<5; i++) for(j=0; j<5; j++)
        {
            (*R)[i][j]=G[1][1][i][j]+G[2][2][i][j];
            (*R)[i][j+5]= -G[1][0][i][j];
            (*R)[i][j+10]= -G[2][0][i][j];
            (*R)[i+5][j]= -G[0][1][i][j];
            (*R)[i+5][j+5]=G[0][0][i][j]+G[2][2][i][j];
            (*R)[i+5][j+10]= -G[2][1][i][j];
            (*R)[i+10][j]= -G[0][2][i][j];
            (*R)[i+10][j+5]= -G[1][2][i][j];
            (*R)[i+10][j+10]=G[0][0][i][j]+G[1][1][i][j];
        }
    return;
}

void local_div_bub_div_bub(const grid3d_t *grid, const elem_t *T, double ***R)
/*
 *  This is the matrix  with entries
 *     R[i][j] = \int (div.b_i)(div.b_j),   0<= i,j < 15,
 *  where
 *     b_i = [phi_i, 0, 0]^T, for  0<= i <5
 *     b_i = [0, phi_i, 0]^T, for  5<= i <10
 *     b_i = [0, 0, phi_i]^T, for 10<= i <15
 */
{
    int i, j;
    double ****G;

    static double B[225], *BB[15]={
        B, B+15, B+30, B+45, B+60,
        B+75, B+90, B+105, B+120, B+135,
        B+150, B+165, B+180, B+195, B+210
    };
    if(R==NULL)
        return;
    else
        *R=BB;
    tensor_grad_bub_grad_bub(grid, T, &G);
    memset(B, 0, sizeof(B));
    for(i=0; i<5; i++) for(j=0; j<5; j++)
        {
            (*R)[i][j]=G[0][0][i][j];
            (*R)[i][j+5]=G[0][1][i][j];
            (*R)[i][j+10]=G[0][2][i][j];
            (*R)[i+5][j]=G[1][0][i][j];
            (*R)[i+5][j+5]=G[1][1][i][j];
            (*R)[i+5][j+10]=G[1][2][i][j];
            (*R)[i+10][j]=G[2][0][i][j];
            (*R)[i+10][j+5]=G[2][1][i][j];
            (*R)[i+10][j+10]=G[2][2][i][j];
        }
    return;
}

void local_stif_bub(const grid3d_t *grid, const elem_t *T, double ***S)
{
    int i, j;
    /*    test tensor_grad_bub_grad_bub
     * static double C[5][5];
     * double ****G;
     */
    static double B[25], *BB[5]={B, B+5, B+10, B+15, B+20};
    double A[4], n[4][3], TV;
    const double W44=(1<<19)/(double) FAC9/3.0;
    if(S==NULL)
        return;
    else
        *S=BB;
    /*    test tensor_grad_bub_grad_bub
     * tensor_grad_bub_grad_bub(grid, T, &G);
     * memset(C,0,sizeof(C));
     * for(i=0;i<5;i++) for(j=0;j<5;j++)
     * {
     * for(int l=0;l<3;l++)
     * C[i][j] += G[l][l][i][j];
     * }
     */
    /* compute area and normal to faces */
    for(i=0; i<4; i++)
        A[i]=ComputeArea(grid->nodes, elem_face_nodes(T->nodes, i), n[i]);
    /* now compute the local stiffness matrix */
    TV=ComputeVolume(grid->nodes, T->nodes);
    for(i=0; i<4; i++)
        for(j=0; j<4; j++)
            (*S)[i][j]=mulvec(n[i], n[j])*A[i]*A[j]/TV/9.0;
    for(i=0; i<4; i++)
        (*S)[i][4]=(*S)[4][i]=0;
    (*S)[4][4]=W44/TV*(sqr(A[0])+sqr(A[1])+sqr(A[2])+sqr(A[3]));
    return;
}

void local_stif_bub_vec(const grid3d_t *grid, const elem_t *T, double ***S)
{
    int i, j;
    static double A[15*15], *AA[15]={
        A, A+15, A+30, A+45, A+60,
        A+75, A+90, A+105, A+120, A+135,
        A+150, A+165, A+180, A+195, A+210
    };
    double **lS;
    if(S==NULL)
        return;
    else
        *S=AA;
    local_stif_bub(grid, T, &lS);
    memset(A, 0, sizeof(A));
    for(i=0; i<5; i++) for(j=0; j<5; j++)
            (*S)[i][j]=(*S)[i+5][j+5]=(*S)[i+10][j+10]=lS[i][j];
    return;
}

void local_bub_grad_bub(const grid3d_t *grid, const elem_t *T, double ***S)
{
    int i, j, l;
    static double B[3][5][5],
            *BB[15]={B[0][0], B[0][1], B[0][2], B[0][3], B[0][4],
        B[1][0], B[1][1], B[1][2], B[1][3], B[1][4],
        B[2][0], B[2][1], B[2][2], B[2][3], B[2][4]};
    double A[4], n[4][3];
    const double W=(1<<9)/(double) FAC7;
    if(S==NULL)
        return;
    else
        *S=BB;
    memset(B, 0, sizeof(B));
    /* compute area and normal to faces */
    for(i=0; i<4; i++)
        A[i]=ComputeArea(grid->nodes, elem_face_nodes(T->nodes, i), n[i]);
    /* (*S)[i   ][j] contains phi_i*d_x phi_j */
    /* (*S)[i+ 5][j] contains phi_i*d_y phi_j */
    /* (*S)[i+10][j] contains phi_i*d_z phi_j */
    for(i=0; i<4; i++)
        for(j=0; j<4; j++)
            for(l=0; l<3; l++)
                B[l][i][j]= -A[j]*n[j][l]/12.0;
    /*{
     *
     * (*S)[i   ][j]=-A[j]*n[j][0]/12.0;
     * (*S)[i+ 5][j]=-A[j]*n[j][1]/12.0;
     * (*S)[i+10][j]=-A[j]*n[j][2]/12.0;
     * }*/
    for(i=0; i<4; i++)
        for(l=0; l<3; l++)
            B[l][4][i]= -(B[l][i][4]=W*A[i]*n[i][l]);
    /*{
     * (*S)[i   ][4]=W*A[i]*n[i][0];
     * (*S)[i+ 5][4]=W*A[i]*n[i][1];
     * (*S)[i+10][4]=W*A[i]*n[i][2];
     * (*S)[4][i]=-(*S)[i   ][4];
     * (*S)[9][i]=-(*S)[i+ 5][4];
     * (*S)[14][i]=-(*S)[i+10][4];
     * }*/
    for(l=0; l<3; l++)
        B[l][4][4]=0.0;
    /*
     * (*S)[4][4]=0.0;
     * (*S)[9][4]=0.0;
     * (*S)[14][4]=0.0;
     */
    return;
}

static void tensor_grad_bub_grad_bub(const grid3d_t *grid, const elem_t *T,
                                     double *****_G)
/*
 *  After the following code is executed
 *     double ****G;
 *     tensor_grad_bub_grad_bub(grid, T, &G);
 *  the returned G is such that
 *
 *                     /
 *     G[a][b][i][j] = |  d_a phi_i . d_b phi_j
 *                     /
 *
 */
{
    static double G[3][3][5][5];
    static double *B[3*3*5]={
        G[0][0][0], G[0][0][1], G[0][0][2], G[0][0][3], G[0][0][4],
        G[0][1][0], G[0][1][1], G[0][1][2], G[0][1][3], G[0][1][4],
        G[0][2][0], G[0][2][1], G[0][2][2], G[0][2][3], G[0][2][4],
        G[1][0][0], G[1][0][1], G[1][0][2], G[1][0][3], G[1][0][4],
        G[1][1][0], G[1][1][1], G[1][1][2], G[1][1][3], G[1][1][4],
        G[1][2][0], G[1][2][1], G[1][2][2], G[1][2][3], G[1][2][4],
        G[2][0][0], G[2][0][1], G[2][0][2], G[2][0][3], G[2][0][4],
        G[2][1][0], G[2][1][1], G[2][1][2], G[2][1][3], G[2][1][4],
        G[2][2][0], G[2][2][1], G[2][2][2], G[2][2][3], G[2][2][4]
    };
    static double **C[3*3]={
        B, B+5, B+10,
        B+15, B+20, B+25,
        B+30, B+35, B+40
    };
    static double ***D[3]={C, C+3, C+6};

    double An[4][3], SG[3], TV, S;
    int a, b, i, j;

    if(_G) *_G=D;
    else return;

    /* compute area and normal to faces */
    for(i=0; i<4; i++)
    {
        S=ComputeArea(grid->nodes, elem_face_nodes(T->nodes, i), An[i]);
        An[i][0]*=S;
        An[i][1]*=S;
        An[i][2]*=S;
    }
    TV=ComputeVolume(grid->nodes, T->nodes);

    SG[0]=SG[1]=SG[2]=0.0;
    for(i=0; i<4; i++) for(a=0; a<3; a++) SG[a]+=An[i][a];

    memset(G, 0, sizeof(G));
    for(a=0; a<3; a++) for(b=0; b<3; b++)
        {
            const double WI4=(1<<8)/(double) FAC5/9.0;
            const double W44=(1<<16)/(double) FAC7/27.0;
            /* case i  (i<4 and j<4) */
            for(i=0; i<4; i++) for(j=0; j<4; j++)
                    G[a][b][i][j]=An[i][a]*An[j][b]/9.0/TV;
            /* case ii (i<4, j=4 or i=4,j<4 */
            for(i=0; i<4; i++)
                G[a][b][i][4]=G[b][a][4][i]=WI4*An[i][a]*SG[b]/TV;
            /* case iii (i=4 and j=4) */
            for(S=i=0; i<4; i++)
                S+=An[i][a]*An[i][b];
            G[a][b][4][4]=W44*(S+SG[a]*SG[b])/TV;
        }

    return;
}

static void tensor_bub_grad_bub_grad_bub(const grid3d_t *grid, const elem_t *T,
                                         double ******_A)
/*
 *  After the following code is executed
 *     double *****H;
 *     tensor_grad_bub_grad_bub(grid, T, &H);
 *  the returned H is such that
 *
 *                        /
 *     H[a][b][i][j][k] = |  phi_i . d_a phi_j . d_b phi_k
 *                        /
 *
 */
{
    static double A[3][3][5][5][5];
    static double *B[3][3][5][5], **C[3][3][5], ***D[3][3], ****E[3];
    int i, j, k, l, a, b;
    double An[4][3], TV;
    const double CONST1=(1<<19)/(double) FAC10;
    const double CONST2=(1<<26)/(double) FAC11/13.0;

    for(a=0; a<3; a++)
    {
        for(b=0; b<3; b++)
        {
            for(i=0; i<5; i++)
            {
                for(j=0; j<5; j++)
                    B[a][b][i][j]=A[a][b][i][j];
                C[a][b][i]=B[a][b][i];
            }
            D[a][b]=C[a][b];
        }
        E[a]=D[a];
    }
    if(_A) (*_A)=E;
    else return;

    /* compute area and normal to faces */
    for(i=0; i<4; i++)
    {
        double S=ComputeArea(grid->nodes, elem_face_nodes(T->nodes, i), An[i]);
        An[i][0]*=S;
        An[i][1]*=S;
        An[i][2]*=S;
    }
    TV=ComputeVolume(grid->nodes, T->nodes);

    memset(A, 0, sizeof(A));
    for(a=0; a<3; a++) for(b=0; b<3; b++)
        {
            double SS=0.0;
            for(l=0; l<4; l++) SS+=An[l][a]*An[l][b];

            /* case j<4, k<4 */
            for(j=0; j<4; j++) for(k=0; k<4; k++)
                {
                    double Val=An[j][a]*An[k][b]/TV;
                    for(i=0; i<4; i++) A[a][b][i][j][k]=Val/36.0;
                    A[a][b][4][j][k]=Val*32.0/945.0;
                }
            /* cases: i<4 with j=4, k<4 or j<4, k=4 */
            for(i=0; i<4; i++) for(l=0; l<4; l++)
                {
                    A[a][b][i][4][l]= -A[a][b][4][i][l];
                    A[a][b][i][l][4]= -A[a][b][4][l][i];
                }
            /* cases: i=4 with j=4, k<4 or j<4, k=4 */
            for(i=0; i<4; i++) for(l=0; l<4; l++)
                {
                    A[a][b][4][4][l]=0.0;
                    A[a][b][4][l][4]=0.0;
                }
            /* case i<4, j=4, k=4 */
            for(i=0; i<4; i++)
                A[a][b][i][4][4]=CONST1*(SS-2.0/3.0*An[i][a]*An[i][b])/TV;
            /* case i=4, j=4, k=4 */
            A[a][b][4][4][4]=CONST2*SS/TV;
        }

}

static void tensor_bub_bub_grad_bub(const grid3d_t *grid, const elem_t *T,
                                    double ****_A)
{
    static double A[3][5][5][5];
    static double *B[3*5*5]={
        A[0][0][0], A[0][0][1], A[0][0][2], A[0][0][3], A[0][0][4],
        A[0][1][0], A[0][1][1], A[0][1][2], A[0][1][3], A[0][1][4],
        A[0][2][0], A[0][2][1], A[0][2][2], A[0][2][3], A[0][2][4],
        A[0][3][0], A[0][3][1], A[0][3][2], A[0][3][3], A[0][3][4],
        A[0][4][0], A[0][4][1], A[0][4][2], A[0][4][3], A[0][4][4],
        A[1][0][0], A[1][0][1], A[1][0][2], A[1][0][3], A[1][0][4],
        A[1][1][0], A[1][1][1], A[1][1][2], A[1][1][3], A[1][1][4],
        A[1][2][0], A[1][2][1], A[1][2][2], A[1][2][3], A[1][2][4],
        A[1][3][0], A[1][3][1], A[1][3][2], A[1][3][3], A[1][3][4],
        A[1][4][0], A[1][4][1], A[1][4][2], A[1][4][3], A[1][4][4],
        A[2][0][0], A[2][0][1], A[2][0][2], A[2][0][3], A[2][0][4],
        A[2][1][0], A[2][1][1], A[2][1][2], A[2][1][3], A[2][1][4],
        A[2][2][0], A[2][2][1], A[2][2][2], A[2][2][3], A[2][2][4],
        A[2][3][0], A[2][3][1], A[2][3][2], A[2][3][3], A[2][3][4],
        A[2][4][0], A[2][4][1], A[2][4][2], A[2][4][3], A[2][4][4]
    };
    static double **C[3*5]={
        B, B+5, B+10, B+15, B+20,
        B+25, B+30, B+35, B+40, B+45,
        B+50, B+55, B+60, B+65, B+70
    };
    double Q[4], n[4][3];
    int i, j, l, k;
    if(_A) *_A=C;
    else return;
    /* compute area and normal to faces */
    for(i=0; i<4; i++)
        Q[i]=ComputeArea(grid->nodes, elem_face_nodes(T->nodes, i), n[i]);
    /* (*_A)[i   ][j][k] contains phi_i**phi_j*d_x phi_k */
    /* (*_A)[i+5 ][j][k] contains phi_i**phi_j*d_y phi_k */
    /* (*_A)[i+10][j][k] contains phi_i**phi_j*d_z phi_k */
    /* case (i): i,j,k < 4 */
    for(i=0; i<4; i++) for(j=0; j<4; j++) for(k=0; k<4; k++)
            {
                for(l=0; l<3; l++)
                    A[l][i][j][k]= -(i==j?4:2)/(double) FAC5*Q[k]*n[k][l];
            }
    /* case (i++) i=4,j<4,k<4 or i<4,j=4,k<4 */
    for(i=0; i<4; i++) for(k=0; k<4; k++)
        {
            const double W=(1<<10)/(double) FAC8;
            for(l=0; l<3; l++)
                A[l][i][4][k]=A[l][4][i][k]= -W*Q[k]*n[k][l];
        }
    for(k=0; k<4; k++)
    {
        const double W=(1<<21)/(double) FAC11;
        for(l=0; l<3; l++)
            A[l][4][4][k]= -W*Q[k]*n[k][l];
    }
    /* case (ii): i,j < 4, k=4 */
    for(i=0; i<4; i++) for(j=0; j<4; j++)
        {
            const double W=(1<<7)/(double) FAC7;
            for(l=0; l<3; l++)
                A[l][i][j][4]=W*(Q[i]*n[i][l]+Q[j]*n[j][l]);
        }
    /* case (iii): i<4, j=4, k=4 or i=4, j<4, k=4 */
    for(i=0; i<4; i++)
    {
        const double W=2*(1<<19)/(double) FAC11;
        for(l=0; l<3; l++)
            A[l][i][4][4]=A[l][4][i][4]=W*Q[i]*n[i][l];
    }
    /* case (iv): i=j=k=4 */
    for(l=0; l<3; l++)
        A[l][4][4][4]=0.0;

    return;
}

void local_advm_bub(const grid3d_t *grid, const elem_t *T, double **AV,
                    double ***_M)
/*
 *  AV must be 3x5, i.e. A = sum_j (A[0][j], A[1][j], A[2][j])^T *phi_j
 */
{
    static double M[5][5];
    static double *MM[5]={M[0], M[1], M[2], M[3], M[4]};
    double ***A;
    int i, j, k, m;

    if(_M) *_M=MM;
    else return;

    tensor_bub_bub_grad_bub(grid, T, &A);
    memset(M, 0, sizeof(M));
    /* This loop fills up G_{ij}^a = int phi_i a.grad phi_j */
    for(i=0; i<5; i++)
        for(j=0; j<5; j++)
            for(m=0; m<3; m++)
                for(k=0; k<5; k++)
                    M[i][j]+=AV[m][k]*A[i+5*m][k][j];
    return;
}

void local_advm_D_bub(const grid3d_t *grid, const elem_t *T, double **AV,
                      double ***_M)
/*
 *  AV must be 3x5, i.e. A = sum_j (A[0][j], A[1][j], A[2][j])^T *phi_j
 */
{
    static double M[5][5];
    static double *MM[5]={M[0], M[1], M[2], M[3], M[4]};
    double ***A;
    int i, j, k, m;

    if(_M) *_M=MM;
    else return;

    tensor_bub_bub_grad_bub(grid, T, &A);
    memset(M, 0, sizeof(M));
    /* This loop fills up D_{ij}^a = 0.5 * int phi_i (div.a) phi_j */
    for(i=0; i<5; i++)
        for(j=0; j<5; j++)
            for(m=0; m<3; m++)
                for(k=0; k<5; k++)
                    M[i][j]+=0.5*AV[m][k]*A[i+5*m][j][k];
    return;
}

void local_advm_bub_vec(const grid3d_t *grid, const elem_t *T, double **AV,
                        double ***M)
{
    int i, j;
    static double A[15][15],
            *AA[15]={A[0], A[1], A[2], A[3], A[4],
        A[5], A[6], A[7], A[8], A[9],
        A[10], A[11], A[12], A[13], A[14]};
    double **lM;
    if(M==NULL)
        return;
    else
        *M=AA;
    memset(A, 0, sizeof(A));
    local_advm_bub(grid, T, AV, &lM);
    for(i=0; i<5; i++) for(j=0; j<5; j++)
            (*M)[i][j]=(*M)[i+5][j+5]=(*M)[i+10][j+10]=lM[i][j];
    return;
}

void local_advm_D_bub_vec(const grid3d_t *grid, const elem_t *T, double **AV,
                          double ***M)
{
    int i, j;
    static double A[15][15],
            *AA[15]={A[0], A[1], A[2], A[3], A[4],
        A[5], A[6], A[7], A[8], A[9],
        A[10], A[11], A[12], A[13], A[14]};
    double **lM;
    if(M==NULL)
        return;
    else
        *M=AA;
    local_advm_D_bub(grid, T, AV, &lM);
    memset(A, 0, sizeof(A));
    for(i=0; i<5; i++) for(j=0; j<5; j++)
            (*M)[i][j]=(*M)[i+5][j+5]=(*M)[i+10][j+10]=lM[i][j];
    return;
}

void local_adjoint_B_bub_vec(const grid3d_t *grid, const elem_t *T, double **AV,
                             double ***_M)
/*
 *  This is the matrix of M[p][q] = int (b_q . grad) AV . b_p ,
 *  where {b_p}_p=0^14 are the vector basis of P1^+
 */
{
    static double M[15][15];
    static double *MM[15]={M[0], M[1], M[2], M[3], M[4],
        M[5], M[6], M[7], M[8], M[9], M[10], M[11], M[12], M[13], M[14]};
    double ***A;
    int i, j, k, a, b;

    if(_M) *_M=MM;
    else return;

    tensor_bub_bub_grad_bub(grid, T, &A);
    /* A[i+5*l][j][k] = int phi_i phi_j d_l phi_k */
    memset(M, 0, sizeof(M));
    // see notes for matrix III with corrections on P_1^+P_1 advection matrices
#if (1)
    for(a=0; a<3; a++) for(i=0; i<5; i++)
            for(b=0; b<3; b++) for(j=0; j<5; j++)
                    for(k=0; k<5; k++)
                        M[i+5*a][j+5*b]+=AV[b][k]*A[i+5*a][j][k];
#else
    for(a=0; a<3; a++) for(i=0; i<5; i++)
            for(b=0; b<3; b++) for(j=0; j<5; j++)
                    for(k=0; k<5; k++)
                        M[i+5*a][j+5*b]+=AV[a][k]*A[j+5*b][i][k];
#endif
}

void local_adjoint_C_bub_vec(const grid3d_t *grid, const elem_t *T, double **AV,
                             double ***_M)
/*
 *  This is the matrix of M[p][q] = int (div b_q) (AV . b_p),
 *  where {b_p}_p=0^14 are the vector basis of P1^+
 */
{
    static double M[15][15];
    static double *MM[15]={M[0], M[1], M[2], M[3], M[4],
        M[5], M[6], M[7], M[8], M[9], M[10], M[11], M[12], M[13], M[14]};
    double ***A;
    int i, j, k, a, b;

    if(_M) *_M=MM;
    else return;

    tensor_bub_bub_grad_bub(grid, T, &A);
    /* A[i+5*l][j][k] = int phi_i phi_j d_l phi_k */
    memset(M, 0, sizeof(M));
    for(a=0; a<3; a++)
        for(i=0; i<5; i++)
            for(b=0; b<3; b++)
                for(j=0; j<5; j++)
                    for(k=0; k<5; k++)
                        M[i+5*a][j+5*b]+=0.5*AV[a][k]*A[i+5*b][k][j];
}

static int __ccc(int i, int j, int k, int *g)
/* returns c and g such that  e_i ^ (e_j ^ e_k) = c e_g */
{
    if(j!=k)
    {
        if(i==j)
        {
            *g=k;
            return -1;
        }
        if(i==k)
        {
            *g=j;
            return 1;
        }
        return 0;
    }
    return 0;
}

void local_adjoint_curl_bub_vec(const grid3d_t *grid, const elem_t *T,
                                double **AV, double ***_M)
/*
 *  This is the matrix of M[p][q] = int (AV x b_q) . (curl b_p),
 *  where {b_p}_p=0^14 are the vector basis of P1^+
 */
{
    static double M[15][15];
    static double *MM[15]={M[0], M[1], M[2], M[3], M[4],
        M[5], M[6], M[7], M[8], M[9], M[10], M[11], M[12], M[13], M[14]};
    double ***A;
    int i, j, k, a, b, c, g;

    if(_M) *_M=MM;
    else return;

    tensor_bub_bub_grad_bub(grid, T, &A);
    /* A[i+5*l][j][k] = int phi_i phi_j d_l phi_k */
    memset(M, 0, sizeof(M));
    for(a=0; a<3; a++) for(b=0; b<3; b++) for(c=0; c<3; c++)
            {
                switch(__ccc(a, c, b, &g))
                {
                    case -1 : /* in this case g == b */
                        for(i=0; i<5; i++) for(j=0; j<5; j++) for(k=0; k<5; k++)
                                    M[i+5*a][j+5*b]-=AV[c][k]*A[j+5*g][k][i];
                        break;
                    case 1: /* in this case g == c */
                        for(i=0; i<5; i++) for(j=0; j<5; j++) for(k=0; k<5; k++)
                                    M[i+5*a][j+5*b]+=AV[c][k]*A[j+5*g][k][i];
                        break;
                    case 0:
                        break;
                    default: error_msg("%s: error in __ccc\n", __FILE__);
                }
            }
}

typedef struct
{
    grid3d_t const *grid;
    field3d *F;
    int el;
    double time;
} func_struct_t;


/****************************************************************************/
/*      Scalar fields in BUB                                                */

/****************************************************************************/

static void func_scal_bub(const double *loc, double **R, void *data)
{
    func_struct_t *d=data;
    int *Tnodes;
    node_t *nodes;
    int i;
    double J, V, B[5], glo[3];
    static double _R[5];

    *R=_R;
    Tnodes=d->grid->elems[d->el].nodes;
    nodes=d->grid->nodes;

    /* the values of the bub basis functions at point loc go in B */
    B[0]=1.0-loc[0]-loc[1]-loc[2];
    B[1]=loc[0];
    B[2]=loc[1];
    B[3]=loc[2];
    B[4]=256*B[0]*B[1]*B[2]*B[3];

    ComputeGlobalCoord(nodes, Tnodes, loc, glo);

    /* the value of the scalar function at the point */
    d->F(d->time, glo, &V);

    /* NOTE: Jacobian  = 6 * Volume */
    J=6.0*ComputeVolume(nodes, Tnodes);
    ;

    for(i=0; i<5; i++)
        _R[i]=J*V*B[i];
}

void test_scalar_field_bub(const grid3d_t *grid, double time, field3d *field,
                           scalar A)
{
    int i;
    double V[5];
    func_struct_t d;
    bub_extend_t *be=get_be(grid);

    if(be->bubNN!=get_N_scal(A))
        error_msg("Wrong size of scalar A in `test_scalar_bub`.\n");

    /* set result to 0 before we begin */
    clean_scal(A);

    /* if field is NULL then the answer is 0 */
    if(field==NULL) return;
    d.grid=grid;
    d.time=time;
    d.F=field;

    use_scalar(A)
    for(d.el=0; d.el<grid->NEL; d.el++)
    {
        integrate_ref_el(QDEG5, 5, func_scal_bub, V, &d);
        for(i=0; i<5; i++)
            AF[be->b_el[d.el].nodes[i]]+=V[i];
    }
    done_scalar(A)
}

void interpolate_scalar_lagr_bub(const grid3d_t *grid, double time, field3d *F,
                                 scalar s)
{
    int i;
    double V;
    assert(s!=NULL&&get_N_scal(s)==ndof_bub(grid));
    clean_scal(s);
    if(F==NULL) return;

    use_scalar(s)
    for(i=0; i<grid->NN; i++)
    {
        F(time, grid->nodes[i].coords, &V);
        sF[i]=V;
    }
    for(i=0; i<grid->NEL; i++)
    {
        elem_t *T=grid->elems+i;
        int j;
        double C[3];
        ComputeCentroid(grid->nodes, T->nodes, C);
        F(time, C, &V);
        for(j=0; j<4; j++)
            V-=0.25*sF[T->nodes[j]];
        sF[i+grid->NN]=V;
    }
    done_scalar(s)
}



/****************************************************************************/
/*      Vector fields in BUB                                                */

/****************************************************************************/

static void func_vect_bub(const double *loc, double **R, void *data)
{
    func_struct_t *d=data;
    node_t *nodes;
    int *Tnodes;
    int i;
    double J, V[3], B[5], glo[3];
    static double _R[15];

    *R=_R;
    nodes=d->grid->nodes;
    Tnodes=d->grid->elems[d->el].nodes;

    /* the values of the bub basis functions at point loc go in B */
    B[0]=1.0-loc[0]-loc[1]-loc[2];
    B[1]=loc[0];
    B[2]=loc[1];
    B[3]=loc[2];
    B[4]=256*B[0]*B[1]*B[2]*B[3];

    ComputeGlobalCoord(nodes, Tnodes, loc, glo);

    /* the value of the vector function at the point */
    d->F(d->time, glo, V);

    /* NOTE: Jacobian  = 6 * Volume */
    J=6.0*ComputeVolume(nodes, Tnodes);
    for(i=0; i<5; i++)
    {
        _R[i ]=J*V[0]*B[i];
        _R[i+5 ]=J*V[1]*B[i];
        _R[i+10]=J*V[2]*B[i];
    }
}

void test_vector_field_bub(const grid3d_t *grid, double time, field3d *field,
                           vector A)
{
    func_struct_t d;
    int i;
    double V[15];
    bub_extend_t *be=get_be(grid);

    if(be->bubNN!=get_N_vec(A))
        error_msg("Wrong size of vector A in `test_vector_bub`.\n");

    /* set result to 0 before we begin */
    clean_vec(A);

    /* if field is NULL then the answer is 0 */
    if(field==NULL)
        return;

    d.F=field;
    d.grid=grid;
    d.time=time;

    use_vector(A)
    for(d.el=0; d.el<grid->NEL; d.el++)
    {
        integrate_ref_el(QDEG5, 15, func_vect_bub, V, &d);
        for(i=0; i<5; i++)
        {
            AX[be->b_el[d.el].nodes[i]]+=V[i];
            AY[be->b_el[d.el].nodes[i]]+=V[i+5];
            AZ[be->b_el[d.el].nodes[i]]+=V[i+10];
        }
    }
    done_vector(A)
}

void interpolate_vector_lagr_bub(const grid3d_t *grid, double time, field3d *F,
                                 vector s)
{
    int i;
    double V[3];

    assert(get_N_vec(s)==ndof_bub(grid));
    clean_vec(s);
    if(F==NULL) return;

    use_vector(s)
    for(i=0; i<grid->NN; i++)
    {
        F(time, grid->nodes[i].coords, V);
        sX[i]=V[0];
        sY[i]=V[1];
        sZ[i]=V[2];
    }
    for(i=0; i<grid->NEL; i++)
    {
        elem_t *T=grid->elems+i;
        int j;
        double C[3];
        ComputeCentroid(grid->nodes, T->nodes, C);
        F(time, C, V);
        for(j=0; j<4; j++)
        {
            V[0]-=0.25*sX[T->nodes[j]];
            V[1]-=0.25*sY[T->nodes[j]];
            V[2]-=0.25*sZ[T->nodes[j]];
        }
        sX[i+grid->NN]=V[0];
        sY[i+grid->NN]=V[1];
        sZ[i+grid->NN]=V[2];
    }
    done_vector(s)
}

void test_grad_bub_in_grad_bub_scal(const grid3d_t *grid, const scalar p,
                                    scalar q)
{
    bub_extend_t *be=get_be(grid);
    double **lD;
    int el;
    assert((q!=NULL)&&(p!=NULL));
    assert((get_N_scal(p)==be->bubNN)&&(get_N_scal(q)==be->bubNN));
    clean_scal(q);

    use_scalar(q)
    use_scalar(p)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, *n;
        local_stif_bub(grid, grid->elems+el, &lD);
        n=be->b_el[el].nodes;
        for(i=0; i<5; i++)
            for(j=0; j<5; j++)
                qF[n[i]]+=lD[i][j]*pF[n[j]];
    }
    done_scalar(p)
    done_scalar(q)
}

void test_grad_bub_in_grad_bub_vec(const grid3d_t *grid, const vector u,
                                   vector v)
{
    bub_extend_t *be=get_be(grid);
    double **lD;
    int el;
    assert((u!=NULL)&&(v!=NULL));
    assert((get_N_vec(u)==be->bubNN)&&(get_N_vec(v)==be->bubNN));
    clean_vec(v);

    use_vector(u)
    use_vector(v)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, *n;
        local_stif_bub(grid, grid->elems+el, &lD);
        n=be->b_el[el].nodes;
        for(i=0; i<5; i++) for(j=0; j<5; j++)
            {
                vX[n[i]]+=lD[i][j]*uX[n[j]];
                vY[n[i]]+=lD[i][j]*uY[n[j]];
                vZ[n[i]]+=lD[i][j]*uZ[n[j]];
            }
    }
    done_vector(v)
    done_vector(u)
}

void test_curl_bub_in_curl_bub_vec(const grid3d_t *grid, const vector u,
                                   vector v)
{
    bub_extend_t *be=get_be(grid);
    double **lD;
    int el;
    assert((u!=NULL)&&(v!=NULL));
    assert((get_N_vec(u)==be->bubNN)&&(get_N_vec(v)==be->bubNN));
    clean_vec(v);
    scalar su=vec2scal(u);
    scalar sv=vec2scal(v);

    use_scalar(su)
    use_scalar(sv)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, *n;
        local_curl_bub_curl_bub(grid, grid->elems+el, &lD);
        n=be->b_el[el].dofs;
        for(i=0; i<15; i++) for(j=0; j<15; j++)
                svF[n[i]]+=lD[i][j]*suF[n[j]];
    }
    done_scalar(sv)
    done_scalar(su)

    free_scal(su);
    free_scal(sv);
}

void test_grad_p1_in_bub(const grid3d_t *grid, const scalar p, vector g)
{
    bub_extend_t *be=get_be(grid);
    double **lD;
    int el;
    assert((g!=NULL)&&(p!=NULL));
    assert((get_N_scal(p)==grid->NN)&&(get_N_vec(g)==be->bubNN));
    clean_vec(g);

    use_scalar(p)
    use_vector(g)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, *n;
        local_bub_grad_bub(grid, grid->elems+el, &lD);
        n=be->b_el[el].nodes; /* the bubble node is always last */
        for(i=0; i<5; i++) for(j=0; j<4; j++)
            {
                gX[n[i]]+=lD[i ][j]*pF[n[j]];
                gY[n[i]]+=lD[i+5][j]*pF[n[j]];
                gZ[n[i]]+=lD[i+10][j]*pF[n[j]];
            }
    }
    done_vector(g)
    done_scalar(p)
}

void test_grad_bub_in_bub(const grid3d_t *grid, const scalar p, vector g)
{
    bub_extend_t *be=get_be(grid);
    double **lD;
    int el;
    assert((g!=NULL)&&(p!=NULL));
    assert((get_N_scal(p)==be->bubNN)&&(get_N_vec(g)==be->bubNN));
    clean_vec(g);

    use_vector(g)
    use_scalar(p)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, *n;
        local_bub_grad_bub(grid, grid->elems+el, &lD);
        n=be->b_el[el].nodes; /* the bubble node is always last */
        for(i=0; i<5; i++) for(j=0; j<5; j++)
            {
                gX[n[i]]+=lD[i ][j]*pF[n[j]];
                gY[n[i]]+=lD[i+5][j]*pF[n[j]];
                gZ[n[i]]+=lD[i+10][j]*pF[n[j]];
            }
    }
    done_scalar(p)
    done_vector(g)
}

void test_p1_in_div_bub(const grid3d_t *grid, const scalar p, vector g)
{
    bub_extend_t *be=get_be(grid);
    double **lD;
    int el;
    assert((g!=NULL)&&(p!=NULL));
    assert((get_N_scal(p)==grid->NN)&&(get_N_vec(g)==be->bubNN));
    clean_vec(g);

    use_vector(g)
    use_scalar(p)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, *n;
        local_bub_grad_bub(grid, grid->elems+el, &lD);
        n=be->b_el[el].nodes; /* the bubble node is always last */
        for(i=0; i<5; i++) for(j=0; j<4; j++)
            {
                gX[n[i]]+=lD[j ][i]*pF[n[j]];
                gY[n[i]]+=lD[j+5][i]*pF[n[j]];
                gZ[n[i]]+=lD[j+10][i]*pF[n[j]];
            }
    }
    done_scalar(p)
    done_vector(g)
}

void test_bub_in_grad_p1(const grid3d_t *grid, const vector u, scalar q)
{

    bub_extend_t *be=get_be(grid);
    double **lD;
    int el;
    assert((u!=NULL)&&(q!=NULL));
    assert((get_N_scal(q)==grid->NN)&&(get_N_vec(u)==be->bubNN));
    clean_scal(q);

    use_vector(u)
    use_scalar(q)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, *n;
        local_bub_grad_bub(grid, grid->elems+el, &lD);
        n=be->b_el[el].nodes; /* the bubble node is always last */
        for(i=0; i<5; i++) for(j=0; j<4; j++)
            {
                qF[n[j]]+=lD[i ][j]*uX[n[i]];
                qF[n[j]]+=lD[i+5][j]*uY[n[i]];
                qF[n[j]]+=lD[i+10][j]*uZ[n[i]];
            }
    }
    done_scalar(q)
    done_vector(u)
}

void test_div_bub_in_p1(const grid3d_t *grid, const vector u, scalar q)
{
    bub_extend_t *be=get_be(grid);
    double **lD;
    int el;
    assert((u!=NULL)&&(q!=NULL));
    assert((get_N_scal(q)==grid->NN)&&(get_N_vec(u)==be->bubNN));
    clean_scal(q);

    use_vector(u)
    use_scalar(q)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, *n;
        local_bub_grad_bub(grid, grid->elems+el, &lD);
        n=be->b_el[el].nodes; /* the bubble node is always last */
        for(i=0; i<5; i++) for(j=0; j<4; j++)
            {
                qF[n[j]]+=lD[j ][i]*uX[n[i]];
                qF[n[j]]+=lD[j+5][i]*uY[n[i]];
                qF[n[j]]+=lD[j+10][i]*uZ[n[i]];
            }
    }
    done_scalar(q)
    done_vector(u)
}

void test_bub_grad_bub_in_grad_p1(const grid3d_t *grid, const vector u,
                                  const vector v, scalar q)
{
    bub_extend_t *be=get_be(grid);
    int el;
    assert(q!=NULL);
    assert(get_N_scal(q)==grid->NN);
    clean_scal(q);
    if(u==NULL||v==NULL) return;
    assert(get_N_vec(u)==be->bubNN&&get_N_vec(v)==be->bubNN);

    use_vector(u)
    use_vector(v)
    use_scalar(q)
    for(el=0; el<grid->NEL; el++)
    {
        int i, j, k, b, c, *dofs=be->b_el[el].nodes;
        double *****lH;
        double U[5][3], V[5][3];
        tensor_bub_grad_bub_grad_bub(grid, grid->elems+el, &lH);
        for(j=0; j<5; j++)
        {
            U[j][0]=uX[dofs[j]];
            V[j][0]=vX[dofs[j]];
            U[j][1]=uY[dofs[j]];
            V[j][1]=vY[dofs[j]];
            U[j][2]=uZ[dofs[j]];
            V[j][2]=vZ[dofs[j]];
        }
        for(b=0; b<3; b++)
            for(c=0; c<3; c++)
                for(j=0; j<5; j++)
                    for(k=0; k<5; k++)
                        for(i=0; i<4; i++)
                            qF[dofs[i]]+=U[j][b]*V[k][c]*lH[b][c][j][k][i];
    }
    done_scalar(q)
    done_vector(v)
    done_vector(u)
}

/****************************************************************************/

void surf_test_bub_dot_n_in_p1(const grid3d_t *grid, const vector u, scalar q)
{
    double A, n[3], Un, U[3];
    int i, j, fid;
    bub_extend_t *be=get_be(grid);
    assert(q!=NULL);
    assert(get_N_scal(q)==grid->NN);
    clean_scal(q);
    if(u==NULL) return;
    assert(get_N_vec(u)==be->bubNN);

    use_vector(u)
    use_scalar(q)
    for(fid=grid->NFI; fid<grid->NF; fid++)
    {
        face_t *f=grid->faces+fid;
        int Tid=f->elems[0];
        elem_t *T=grid->elems+Tid;
        int s=0;
        for(s=0; s<4; s++) if(T->faces[s]==fid) break;
        A=ComputeArea(grid->nodes, elem_face_nodes(T->nodes, s), n);
        for(j=0; j<3; j++)
        {
            U[0]=uX[f->nodes[j]];
            U[1]=uY[f->nodes[j]];
            U[2]=uZ[f->nodes[j]];
            Un=mulvec(U, n);
            for(i=0; i<3; i++)
                qF[f->nodes[i]]+=Un*A/(i==j?6.0:12.0);
        }
    }
    done_scalar(q)
    done_vector(u)
}

void surf_test_n_cross_curl_bub_in_grad_p1(const grid3d_t *grid,
                                           const vector u, scalar q)
{
    double TV, A[4], n[4][3], U[5][3];
    int i, j, b, s, fid, Tid, *dofs;
    bub_extend_t *be=get_be(grid);
    assert(q!=NULL);
    assert(get_N_scal(q)==grid->NN);
    clean_scal(q);
    if(u==NULL) return;
    assert(get_N_vec(u)==be->bubNN);

    use_scalar(q)
    use_vector(u)
    for(fid=grid->NFI; fid<grid->NF; fid++)
    {
        face_t *f=grid->faces+fid;
        Tid=f->elems[0];
        elem_t *T=grid->elems+Tid;
        dofs=be->b_el[Tid].nodes;
        for(i=0; i<4; i++)
            A[i]=ComputeArea(grid->nodes, elem_face_nodes(T->nodes, i), n[i]);
        TV=ComputeVolume(grid->nodes, T->nodes);
        for(i=0; i<5; i++)
        {
            U[i][0]=uX[dofs[i]];
            U[i][1]=uY[dofs[i]];
            U[i][2]=uZ[dofs[i]];
        }
        for(s=0; s<4; s++) if(T->faces[s]==fid) break;
        if(s>=4) error_msg("%s: this is impossible.\n", __func__);
        for(i=0; i<4; i++)
        {
            double nis[3], njb[3], eb[3];
            if(i==s) continue;
            crossp(n[i], n[s], nis);
            for(b=0; b<3; b++)
            {
                eb[0]=eb[1]=eb[2]=0.0;
                eb[b]=1.0;
                for(j=0; j<4; j++)
                {
                    crossp(n[j], eb, njb);
                    qF[T->nodes[i]]+=
                            U[j][b]*A[i]*A[j]*A[s]/9.0/TV/TV*mulvec(njb, nis);
                }
                crossp(n[s], eb, njb);
                qF[T->nodes[i]]+=
                        U[4][b]*A[i]*A[s]*A[s]*64.0/135.0/TV/TV*mulvec(njb, nis);
            }
        }
    }
    done_vector(u)
    done_scalar(q)
}


